Thermalization of an anisotropic granular particle 
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We investigate the dynamics of a needle in a two-dimensional bath composed of thermalized point 
particles. Collisions between the needle and points are inelastic and characterized by a normal resti- 
tution coefficient a < 1. By using the Enskog-Boltzmann equation, we obtain analytical expressions 
for the translational and rotational granular temperatures of the needle and show that these are, 
in general, different from the bath temperature. The translational temperature always exceeds the 
rotational one, though the difference decreases with increasing moment of inertia. The predictions 
of the theory are in very good agreement with numerical simulations of the model. 
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I. INTRODUCTION 

The dissipative nature of the collisions in granular sys- 
tems leads to fundamentally different behavior from their 
thermal analogues. A striking example is the lack of 
energy equipartition between the degrees of freedom in 
granular systems. For example experiments 0, Q] and 
computer simulations FJ. 01 have shown that in binary 
mixtures of isotropic inelastic particles, the granular tem- 
peratures of the two species are different. 

These results have prompted theoreticians to investi- 
gate some simple model systems. For example, Mar- 
tin and Piasecki 6] examined the behavior of a spheri- 
cal tracer particle immersed in a homogeneous fluid in 
equilibrium at a temperature T. They showed that the 
Enskog-Boltzmann equation of the tracer particle pos- 
sesses a stationary Maxwellian velocity distribution char- 
acterized by an effective temperature that is smaller than 
T. 

Although many granular systems contain particles that 
are manifestly anisotropic, most studies have been con- 
fined to spherical particles. A notable exception is the 
work of Aspelmeier, Huthmann and Zippelius 0,01 that 
examines the free cooling of an assembly of inelastic nee- 
dles with the aid of a pseudo-Liouville operator. They 
predicted an exponentially fast cooling followed by a state 
with a stationary ratio of translational and rotational 
energy. This two stage cooling was confirmed by event 
driven simulations. 

In this work we examine the breakdown of equipar- 
tition in a steady state granular system containing an 
anisotropic particle. Motivated by the studies of Mar- 
tin and Piasecki[(^ and Aspelmeier, Huthmann and 
Zippelius0, 0, we consider an infinitely thin inelastic 
needle immersed in a bath of point particles. Start- 
ing from the collision rules of this model, we derive 
the Enskog-Boltzmann equation. By assuming that the 
points are thermalized and that the velocity and angu- 
lar velocity distributions of the needle are Maxwellian, 
we derive analytical expressions for the translational and 
rotational granular temperatures as a function of the 



masses of the two species, the moment of inertia of the 
needle and the normal restitution coefficient. Both these 
temperatures are smaller than that of the bath. The ro- 
tational granular temperature is usually lower than the 
translational one, except for very light bath particles for 
which the two are equal. 

We also report essentially exact numerical results ob- 
tained using a stochastic simulation method. The theory 
is in very good agreement with the simulation, which 
validates the assumption of the Maxwellian shape of the 
distribution functions. 



II. MODEL 

We examine a two-dimensional system consisting of an 
infinitely thin needle of length L, mass M and moment 
of inertia I immersed in a bath of point particles each 
of mass m. The vector position of the center of mass of 
the needle and a point particle are denoted by ri and r2 , 
respectively. The orientation of the needle is specified by 
a unit vector Ui that points along its axis. Let ri2 = 
ri — Y2 and denote a vector perpendicular to Ui. A 
collision between a needle and a point occurs when 



ri2.ii! = 



(1) 



and |A| < L/2. (see Fig. QJ. The relative velocity of the 
point of contact V is given by 



V = V12 + Aiii. 



(2) 



The pre- and post-collisional quantities (the latter are 
labeled with a prime) obey the usual conservation laws: 



• Total momentum conservation 

Mv'i + TOv'2 = Afvi + mv2. 



(3) 



• Angular momentum conservation with respect to 
the point of contact 



Iu[k = Iuik + MAui A (v'i - vi) 



(4) 
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FIG. 1: Geometry of the needle and a point in the plane: ri2 
denotes a vector joining the point labeled 2 and the center of 
the needle; ui is a unit vector along the axis of the needle, A 
is the algebraic distance between the center of the needle and 
the point of impact and Ui" is a unit vector perpendicular to 
the axis of the needle. For a collision to occur one requires 
that |A| < L/2 when the point lies on the line defined by the 
needle; i.e. when Eq.Q is satisfied. 



where k is a unit vector perpendicular to the plane 
such that k = ui Au|. 



As a result of the collision, the relative velocity of the 
contacting points changes instantaneously according to 
the following relations: 

V'.u^ = -aV.u| (5) 
V'.ui=V.ui (6) 

where a denotes the normal restitution coefficient. For 
the sake of simplicity we have taken the tangential resti- 
tution coefficient equal to one. This choice is reflected in 
the form of Eq. . 

By combining Eqs. J2J- © one obtains, after some 
algebra, the change of the needle momentum Ap = 
M(vi- Vl ) 

m + M / 

The inelastic collision leads to a loss of translational ki- 



III. PSEUDO-LIOUVILLE EQUATION 

For particles with hard-core interactions, the kinetic 
evolution of iV-particle distribution f(r N ,v N ) is de- 
scribed by a pseudo-Liouville operator, where is a 
short-hand notation for the positions (and internal de- 
grees of freedom, i.e. angle 6i with the s-axis for the 
needle) of N particles and for their velocities (and 
angular velocities for the needle). Originally derived 
for spheres and elastic collisions the formalism has 
been extended to systems composed of inelastic particles. 
Since collisions are instantaneous [Tol Hl| . the pseudo- 
Liouville operator L(r N ,v N ) = L°(r N ) + J2 l<J T lj is the 
sum of the free particle streaming operator L°(r N ) and 
of the sum of two-body collision operators Tij . Since we 
are interested in the homogeneous state, the distribution 
function /(vi, uj\) of the needle system is described by 

-Q t =N / — / dv 2 / dr 2 Ti 2 /(vi, u u v 2 ) 

(10) 

where N is the total number of points, /(vi,wi,v 2 ) the 
distribution function of the needle and a point, and Ti 2 
is the collision operator between a needle and a point. 

The assumption of molecular chaos yields a factor- 
ization of the two body distribution /(vi, u>x, v 2 ) = 
/(vi, lji)$(v 2 ), where $(v 2 ) denotes the point distribu- 
tion function. Since the point particles are thermalized, a 
stationary state for the system (tracer needle and points) 
can be reached, and one finally obtains for the needle dis- 
tribution 

J d6x J dv 2 J dr a r 12 /(vi,wi)$(va) =0 (11) 
where 

$(v 2 ) ex exp (12) 
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where T is the temperature of the bath. 

To build the point-needle collision operator, T12, one 
must include the change in quantities (i.e. velocity and 
angular momentum) produced during the infinitesimal 
time interval of the collision. This operator is different 
from zero only if the two particles are in contact and if 
the particles were approaching just before the collisionQ. 
The explicit form of the operator is 



T 12 cxe(L/2- |A|)5(|r 12 .u^| - 0+ 
d\ri2-ui\ ( d\r 12 .uf- 



dt 



dt 



(&12-I), (13) 



where 612 is an operator that changes pre-collisional 
quantities to post-collisional quantities and Q(x) is the 
Heaviside function. 

The other terms of the collision operator corre- 
spond to the necessary conditions of contact 0(L/2 — 

|A|)£(|r 12 .if I - 0+), and approach 6 (- | d|ri ^' |). 

For an isotropic tracer particle in an isotropic par- 
ticle bath, Martin and PiaseckT6| solved the station- 
ary Enskog-Boltzmann equation, showing that the ve- 
locity distribution of the tracer particle remains gaus- 
sian in a thermalized bath. It is worth noting that 
for finite dilutions, the velocity distribution function is 
not MaxwellianQ, but the deviations from the gaussian 
shape of the distribution function of the velocities (that 
can be calculated in a perturbative way) yields small 
corrections to the estimate of the granular temperature. 
More significant deviations are expected for finite dilu- 
tion systems [l2^ . 

For a mixture of a needle and points, one can show 
that a Maxwell distribution of the needle angular and 
translational momenta cannot be a solution of the sta- 
tionary Enskog-Boltzmann equation even for the case of 
an infinitely diluted needle. This is due to the fact that 
the change of momentum depends on the location of the 
point of impact on the needle (the right-hand side of 
Eq. (JJJ depends on A). However, we impose gaussian 
distributions for the translational and angular velocities 
of the needle. We show below that this is a very good 
approximation. 

In the stationary state, the loss of translational and 
rotational kinetic energies is on average zero and can be 
expressed by means of the collision operator 



(AEf) = (Ti 2 /(vi,wi)$(v 2 )£f (vi)) = 

(ABf) = (Ti 2 /(vi,^i)$(v 2 )i;f (wi)) = (14) 

where the brackets denote an average over the indepen- 
dent variables. This corresponds to taking the second 
moments of the distribution function of the needle. Af- 
ter substitution of the collision operator (Eq.JT^J)), one 
obtains explicitly for the translational kinetic energy 



J ...J dr 2 dv 1 dv 2 cL) 1 d0 1 Q(L/2-\\\)5(\r 1 2.ui\-O + ) 



d|r 12 .u^-| 



dt 



e - 



d\r 12 .u^\ 



dt 



/(v 1 ,^ 1 )$(v 2 )A£; l T = 0. 

(15) 



A similar equation can be written for rotational kinetic 
energy. The needle distribution function f(yi,u>x) is then 
given by 



/(vi,wi) oc exp 



Mvi 5 



(16) 



where 7t and 7^ are the ratios of the translational and 
rotational needle temperatures to the bath temperature, 
respectively. 

By inserting Eq. (|16fl in Eq. I|15fl and in the correspond- 
ing equation for the rotational energy one obtains, after 
some tedious but straightforward calculation (outlined in 
Appendix A), the following set of equations 



1 \/l + akx 2 
axb- 



1 



1 + kx 2 



1 2 \fl + akx 2 

axax ; — 3 — 

1 + kx 2 



1 (1 + akx 2 ) 3 / 2 
o (1 + kx 2 ) 2 ' 



dx 



where 



41 (- + — V 



(17) 

x 2 (l + akx 2 f/ 2 
(1 + kx 2 ) 2 ' 
(18) 



(19) 



and 



a = Jr 



b = 7t 



(M + m) 
M + m^T ' 



M + m 
M + mrjT 



(20) 



(21) 



Explicit expressions for the integrals appearing in 
Eqs. ljT7J) and (fT5f) are given in Appendix B. Eq. 
is an implicit equation for a that, for a given value of 
a, can be solved with standard numerical methods, b is 
then easily obtained by calculating the ratio of integrals 
of Eq. (|17fl . Finally, from the values of a and b, 7t and 
7_r can be obtained from Eas. lf2U|l - !j?T|) . 



IV. RESULTS 

We first note that for the elastic case, a = 1, we have 
the solution a = b = 1 which gives 7t = 1 and since 
a/b = 7r/7t . 7 p = 1 which corresponds to the equilib- 
rium case [l3lll4| . 
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FIG. 2: Ratio of the translational (full curve) 7t and rota- 
tional 7_r (dashed curve) granular temperature to the tem- 
perature of the bath versus the normal restitution coefficient 
a for a homogeneous needle with M = m. 

Figure [5] shows the ratio of the translational tempera- 
ture to the bath temperature and of the rotational tem- 
perature to the bath temperature for a homogeneous nee- 
dle (7 = ML 2 /12) whose mass is equal to the mass of a 
bath particle. As the normal restitution coefficient de- 
creases from 1 to both ratios decrease monotonically 
from 1 to a strictly positive value, such that the trans- 
lational temperature is always larger the rotational tem- 
perature. An analogous situation has been noted for the 
free cooling of needles in three dimensions 7] . 

Figure |3| shows the same ratios for an inhomogeneous 
needle (I = ML 2 /16) corresponding to relatively lighter 
ends and a heavier center. Again the translational tem- 
perature is always larger than the rotational tempera- 
ture, and the differences are greater than for a homo- 
geneous needle. The difference between the two tem- 
peratures vanishes in the limit of very large moment of 
inertia. Contrary to the free cooling state of needles Q, 
there is no "critical" value of the moment of inertia above 
which the rotational temperature is larger than the trans- 
lational temperature. 

We now return to consideration of a homogeneous nee- 
dle and investigate the effect of varying the mass ratio, 
m/M, at constant normal restitution coefficient. When 
m/M — > 0, k — ► 0, which gives a = b = (l + a)/2. Hence, 
for very light bath particles, equipartition between trans- 
lational and rotational granular temperatures is asymp- 
totically obtained and we have 

1 + a 

7fi = 7t = ^— ■ (22) 
This same limit is obtained for a spherical tracer particle 




a 



FIG. 3: Ratio of the translational (full curve) 7t and rota- 
tional 7_r (dashed curve) granular temperature to the tem- 
perature of the bath versus the normal restitution coefficient 
a for an inhomogeneous needle (7 = ML 2 /I6)with M = m. 

. We conjecture that the result is independent of the 
shape of the tracer particle and the dimension. 

Figure 0] displays the variation with the mass ratio of 
the translational and rotational temperature of the tracer 
needle. Note that when the mass of the bath particles 
is very small compared to the needle particle, the two 
curves go to the limit given by Eq. 1221) . 

V. SIMULATION 

In order to assess the accuracy of the theory we de- 
veloped a Direct Simulation Monte Carlo (DSMC) code 
to obtain numerical solutions of the non-linear homoge- 
neous Boltzmann equation. The simulation generates a 
sequence of collisions with a variable time interval be- 
tween each. 

At the beginning of each step, we first calculate the 
total flux of colliding particles on both sides of the nee- 
dle which is a function of the component of the velocity 
of the center of mass normal to the length of the nee- 
dle Vi-uj^ and the angular velocity ui. Next a waiting 
time to the next collision is sampled from an exponential 
distribution with rate constant equal to the total flux. 
The needle is rotated at its current angular velocity to 
the point of collision. The location of the collision is 
then selected with a probability proportional to the po- 
sition dependent flux. The ratio of the flux on the left 
hand side to the total flux is then compared to a uni- 
form random number between zero and one. If the ratio 
is greater than the random number, the collision occurs 
on the right; otherwise is occurs on the left. The normal 



6] 
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FIG. 4: Ratio of the translational (full curve) 7t and rota- 
tional 7_r (dashed curve) granular temperature to the tem- 
perature of the bath versus the ratio of the masses m/M for 
a homogeneous needle (7 = ML 2 / 12 ) and a = 0.9. 

component of the colliding point is then selected from a 
gaussian flux. Finally, the collision rules are applied and 
the normal velocity and angular velocity of the needle arc 
updated. 

The simulation results for 10 7 collisions per run are 
compared with the theory in Figs. |3 Although the 
theory is not exact, the agreement is excellent. The dis- 
crepancy between simulation and theory is never greater 
than 1% whatever the value of the restitution coefficient. 
In Fig. 01 when the mass of the needle increases, the gran- 
ular temperatures obtained by simulation are slightly 
smaller than the analytical result. 

Fig. [5] and Fig. [5] show the translational and angular 
velocity distributions, respectively, for two different val- 
ues of the normal restitution coefficient (0.5 and 0.9). 
On the scale of the plot, the simulated distributions are 
in quantitative agreement with gaussian distributions at 
the corresponding theoretical granular temperature, con- 
firming the validity of employing the latter in the theory. 

VI. CONCLUSION 

We have shown that the translational and rotational 
granular temperatures of an anisotropic tracer particle 
immersed in a bath of point particles depend on the ra- 
tio of the masses of a point and the needle and the mo- 
ment of inertia. They are, in general, not equal and dif- 
fer from the bath temperature. In addition, we found 
that equipartition is obtained asymptotically, regardless 
of the normal restitution coefficient, for very light bath 
particles. We expect this to be a general feature, i.e. in 



FIG. 5: Distribution of the normal velocity component of 
the needle for two values of the normal restitution coefficient, 
a = 0.5 and a = 0.9 (broader curve). The solid (noisy) curves 
are the simulation results and the dashed curves correspond 
to Gaussian distributions of velocity evaluated for the corre- 
sponding theoretical granular temperature. 




FIG. 6: Same asQ3 except the angular velocity is shown. 



higher dimensions and for arbitrary shaped anisotropic 
particles. It should be possible to test this prediction 
experimentally. 
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APPENDIX A: NEEDLE AVERAGE ENERGY 
LOSS 

As for binary mixtures of spheres Q, we use a gaussian 
ansatz for the distribution functions and introduce two 
different temperatures corresponding to the translational 
and rotational degrees of freedom of the needle. The 
homogeneous distribution functions of the needle and of 
the points are then given respectively by 



/(vi,wi) ~ exp 



Mv?7 T 1 Iu{ lR 



2„,-l 



2T 



2T 



<fr(v 2 ) ~ exp 



2T 



(Al) 



(A2) 



where T is the temperature of the bath, 7^ and 7^ the 
ratio of the translational (and rotational) temperature of 
the needle to the bath temperature. 

We introduce the vectors x an d v such that 



1 



V / 2T(M 7T + to) 
/ toM 



(Mvi + TOV 2 ) 



(vi + 7 T v 2 ) 



y 2T{M^(t + to)7t 
The scalar product V.u;^ can be expressed as 



(A3) 
(A4) 



V.uj 1 = 



2T 



[(7t - l)X-ui 



IM 



HT \\hof + V — ^ u i 
1 " M V to / 



(A5) 



Let us introduce £ = ^1 y 2T7r ' ^ e translational energy 
loss is given by the formula 



^ J d\ J d6 x J d X J dv J duj x 



p=±i 



L 



exp (-x 2 - - £ 2 ) |V.u^|e(pV.u^)e (JA| - - J AEf 

(A6) 

Since Eq. (|A5(I depends only on X- U i~ an d is.uj- , one can 
freely integrate over the direction of Ui for the vectors x 
and v. The integration over Q\ can be easily performed. 
If we introduce the three dimensional vectors G and s 
with components: 



G — (Gi, G 2 , G3) 



2T 



Mjt + rn 



( 7t - 1), 



' 2T 7T to + M 
Mj T + m y/mM ' 



(A7) 



and 



s = (81,82,53) 

(x-ui,v.u^,£) 



(A8) 



By inserting Eq. (JBJ in Eq. (|A6|I , the average energy loss 
can be rewritten as 

J d\ J dsexp (-s 2 ) |G.s|6(pG.s)6 ^|A| - | 



2T 



(1 + a)G.s 

(l + a) 2 (G.s) 2 
_ 2M(i + i + ^) 



7TS1 



TO7T 



-82 



(A9) 



By defining a new coordinate system in which the z-axis 
is parallel to G, one find that the integrals of Eq. (|A9|) 
involve gaussian integrals of the form 

j dsexp (-s 2 ) (|G| Sz ) 2 e(±s,)G iS , = ^|G| 2 G, (A10) 



and 



j dsexp (-s 2 ) (|G| S2 ) 3 8(± S2 ) = ^|G 3 | (All) 

which finally leads to Eq. (|17() . The equation for rota- 
tional energy is derived following exactly the same pro- 
cedure. 



APPENDIX B: INTEGRALS 

The integrals appearing in Eqs. (|17J) and Ijl8(l can be 
evaluated explicitly. For completeness, we give below 
their analytical expressions 



h 



VT + akx 2 

dx— , „ 

1 + kx 2 



- ln( V ak + VI + ak) 



1 - a 



■ arctan 



(1 - a)k 
l + ak 



The integral I 2 is defined as 

* Jl +<**>)*'* 



dx- 



, (1 + kx 2 ) 2 
and satisfies the equation 



h=h-{a- l)k 



dd 



(Bl) 



(B2) 



(B3) 



d=k 



7 



which gives 



Finally, one has 



h = — 1= MS ok + VI + ok) + - — n ' 

Vfc 2(1 + fc) 



1 + a - 2a 1 
2y/k{l - a) 



i arctan 



(1 - a)k 
l + ak 



(B4) 



dx- 



2 Vl + akx 2 



1 + fcx 2 



dx\/ 1 + akx 2 — I] 



Vafc+T 1 - 2a 
2fc vafc 3/2 



+ ln(V ak + vT+afc) 



f.3/2 



■ arctan 



(1 - a)k 
l + ak 



(B5) 



(B6) 



Ia = I dx 



1 o; 2 (l + a/fcx 2 ) 3 / 2 



/„ (1 + fcx 2 ) 2 

which satisfy the algebraic equation 

7 2 + = aJi + {a - l)Ii 

which gives 



/4 



VI + afc(afe - 1 + 2a) 
2k{k + 1) 



Vl^a(l - 4a) 
2137 2 



arctan 



(1 - a)fc 
1 + aA: 



(B7) 



(B8) 



(B9) 
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